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Abstract 

We  propose  a  novel  method  for  alleviating  the  stringent  CFL  condition  imposed  by 
the  sound  speed  in  simulating  inviscid  compressible  flow  with  shocks,  contacts  and 
rarefactions.  Our  method  is  based  on  the  pressure  evolution  equation,  so  it  works 
for  arbitrary  equations  of  state,  chemical  species  etc,  and  is  derived  in  a  straight¬ 
forward  manner.  Similar  methods  have  been  proposed  in  the  literature,  but  the 
equations  they  are  based  on  and  the  details  of  the  methods  differ  significantly. 
Notably  our  method  leads  to  a  standard  Poisson  equation  similar  to  what  one 
would  solve  for  incompressible  flow,  but  has  an  identity  term  more  similar  to  a 
diffusion  equation.  In  the  limit  as  the  sound  speed  goes  to  infinity,  one  obtains 
the  Poisson  equation  for  incompressible  flow.  This  makes  the  method  suitable  for 
two-way  coupling  between  compressible  and  incompressible  flows  and  fully  implicit 
solid-fluid  coupling,  although  both  of  these  applications  are  left  to  future  work.  We 
present  a  number  of  examples  to  illustrate  the  quality  and  behavior  of  the  method 
in  both  one  and  two  spatial  dimensions,  and  show  that  for  a  low  Mach  number  test 
case  we  can  use  a  CFL  number  of  300  (whereas  previous  work  was  only  able  to  use 
a  CFL  number  of  3  on  the  same  example). 


1  Introduction 


In  this  paper,  we  focus  on  highly  nonlinear  compressible  flows  with  shocks, 
contacts  and  rarefactions,  for  example  the  Sod  shock  tube.  Traditionally  these 
types  of  problems  are  solved  with  explicit  time  integration  (for  example  Runga- 
Kutta  methods,  ENO,  WENO  etc,  see  e.g.  [10,11,5]).  Although  these  methods 
produce  high  quality  results,  small  time  steps  are  required  in  order  to  enforce 
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the  CFL  condition  of  information  moving  only  one  grid  cell  per  time  step. 
While  this  is  understandable  for  very  high  Mach  number  flow  where  |m|,  Im  — c| 
and  \u  +  c|  are  all  of  similar  magnitude,  it  is  too  restrictive  for  flows  where 
c  may  be  much  larger  than  |m|.  Moreover  some  flow  helds  might  have  both 
high  Mach  number  regions  where  shock  waves  are  of  interest  as  well  as  low 
Mach  number  regions  where  the  material  velocities  are  important.  In  this  case, 
a  large  number  of  time  steps  are  required  if  one  is  interested  in  the  motion 
of  the  fluid  particles  over  an  appreciable  distance  in  the  low  Mach  number 
regions.  Thus,  it  can  be  quite  useful  to  have  methods  that  avoid  the  stringent 
CFL  time  step  restriction  imposed  by  the  acoustic  waves  and  instead  use  only 
the  material  velocity  CFL  restriction  (albeit  one  would  expect  some  loss  of 
quality  because  of  the  implicit  treatment  of  the  acoustic  waves). 


To  alleviate  the  stringent  CFL  restriction,  [6]  proposed  both  a  non-conservative 
and  a  conservative  scheme.  Their  non-conservative  scheme  builds  on  the  predictor- 
corrector  type  scheme  of  [16]  to  derive  an  elliptic  pressure  equation  quite  sim¬ 
ilar  to  ours,  but  for  an  adiabatic  fluid.  Our  method  is  similar  in  spirit  to 
[6,13-15]  where  we  divide  the  calculation  into  two  parts:  advection  and  non- 
advection.  The  advection  terms  are  treated  with  explicit  time  integration,  and 
thus  the  CFL  restriction  on  the  material  velocity  remains.  Whereas  one  can 
use  a  standard  method  such  as  ENO  in  solving  the  advection  terms,  we  found 
that  when  coupled  to  an  implicit  solution  of  the  pressure  equations  (that  is 
inherently  central-differenced)  the  standard  ENO  method  sometimes  leads  to 
spurious  oscillatory  behavior.  Thus  we  designed  a  new  ENO  method  geared 
towards  a  MAC  grid  discretization  of  the  data,  making  it  more  similar  to 
incompressible  flow.  We  call  this  MAC-ENO  or  MENO.  The  remaining  non- 
advection  terms  are  solved  using  an  implicit  equation  for  the  pressure  using 
a  standard  MAC  grid  type  formulation.  Since  the  MAC  grid  is  dual  in  both 
velocity  and  pressure  (noting  that  the  MAC  grid  pressure  needs  to  live  at  cell 
faces  for  flux  based  methods),  one  needs  to  interpolate  data  back  and  forth. 


We  base  the  derivation  of  our  method  on  the  pressure  evolution  equation  as 
discussed  in  [2],  thus  making  it  valid  for  general  equations  of  state,  arbitrary 
chemical  species  etc.  Thus,  our  derivation  has  less  assumptions  and  is  a  bit 
more  straight  forward  than  previous  work,  especially  that  based  on  precon¬ 
ditioners.  Also,  our  method  is  fully  conservative  and  thus  shocks  are  tracked 
at  the  right  speed.  We  present  a  number  of  traditional  examples  for  highly 
non-linear  compressible  flows  including  the  Sod  shock  tube,  interacting  blast 
waves,  and  two  dimensional  flow  past  a  step.  We  also  demonstrate  that  the 
method  works  well  for  low  Mach  number  flow,  taking  the  example  of  [7]  where 
the  authors  obtain  reasonable  results  with  a  CFL  number  of  3.  Notably,  our 
method  allows  a  CFL  number  of  300  (two  orders  of  magnitude  more). 
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2  Numerical  Method 


Let  us  consider  the  one  dimensional  Euler  equations, 
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with  p  being  the  density,  u  the  velocity,  E  the  total  energy  and  p  the  pressure. 
The  flux  term  can  be  separated  into  an  advection  part  and  a  non-advection 
part, 
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We  first  compute  the  advection  part  with  Jacobian 
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All  the  Jacobian’s  eigenvalues  are  equal  to  u,  and  it  is  rank  deficient  with  left 
eigenvectors  of  (u,  — 1,0)  and  (E/p,  0,-1)  and  right  eigenvectors  of  (1,m,0)^ 
and  (0,  0, 1)^.  Since  all  the  characteristic  velocities  are  identical,  we  can  apply 
component  wise  upwinding  to  Fi(U)  without  having  to  transform  into  the 
characteristic  variables  first  (as  in  [4]).  Moreover,  this  advection  part  only 
requires  a  time  step  restriction  based  on  u. 


2.1  Implicit  Pressure  Update 


The  multi-dimensional  Euler  equations  are 
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,  where  u  =  {u,v,w)  are  the  velocities.  Here  we  have  advection  components 
in  each  of  the  3  spatial  dimensions,  and  they  can  be  handled  as  outlined 
previously  in  a  dimension  by  dimension  fashion  as  in  [11], 


We  apply  a  time  splitting  as  is  typical  in  incompressible  flow  formulations,  hrst 
updating  the  advection  terms  to  obtain  an  intermediate  value  of  the  conserved 
variables  (p)*,  {pu)*,  and  E*,  and  afterward  correct  these  to  time  using  the 
pressure.  Since  the  pressure  does  not  affect  the  continuity  equation  =  p*. 
The  non  advection  momentum  and  energy  updates  are 


{pu)^+^  -  {pu)* 

At 


— Vp 
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and 
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At 


=  —V  ■  (pu). 
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As  in  [2],  we  can  use  the  Euler  equations  to  derive  the  pressure  evolution 
equation, 

Pt  +  u  ■  Vp  =  —pc^V  ■  u.  (4) 


Taking  motivation  from  standard  incompressible  flow  solvers,  we  £x  V  ■  u  to 
be  at  time  n  +  1  through  the  time  step,  making  an  0{At)  error. 


Pt  +  u-  Vp  =  —pc^V  ■ 


(5) 


Dividing  equation  (2)  by  p^^^,  and  noting  that  p*  =  p”^^,  gives 

^+1  =  e  -  At^  (6) 

pn+l  ^  ^ 

,  and  following  a  typical  derivation  of  incompressible  flow  we  take  the  diver¬ 
gence  of  equation  (6)  to  obtain 
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Substituting  this  into  equation  (5)  gives 

pt  +  u-Vp=  -pc^V  ■  H*  +  pc^AtV  ■  >  (8) 

which  is  an  advection-diffusion  equation  with  source  term.  Discretizing  the  u  ■ 
Vp  advection  term  explicitly,  using  a  forward  Euler  time  step,  and  defining  the 
diffusive  pressure  at  time  as  is  typical  for  backward  Euler  discretization, 
gives  after  rearrangement 

pH+l  _  =  {p^  -  (PT  ■  Vp”)  At)  -  p^^C^AtV  ■  it.  (9) 


4 


Note  we  have  discretized  p(?  at  time  t".  This  equation  can  be  further  simplihed 
by  using  the  advection  equation  for  pressure, 


At 


+  ir-vp^ 
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to  obtain 

.  \/p^)At, 

where  is  an  advected  pressure  which  can  be  computed  using  HJ  ENO  [9] 
or  semi-Lagrangian  advection  [1],  Substituting  in  equation  (9)  we  obtain 
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We  discretize  this  equation  at  cell  centers,  which  is  typical  for  advection- 
diffusion  equations,  and  thus  need  to  dehne  velocities  at  cell  faces  for  V  ■  n*. 
Consider  two  adjacent  grid  cells  one  centered  at  Xj  and  one  centered  at  Xj+i. 


Fig.  1. 

We  divide  these  into  four  regions  Q^r,  Q+ip,  Ci+i^R,  where 
represents  a  dual  cell  (see  hgure  1).  Then  equation  (2)  for  Ci^R  is 
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Similarly  for  Cj+ip  we  have 
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Adding  these  equations  together  and  dividing  by  (p*  +  Pi+i)  yields 
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where  Ui+1/2  =  pg  thought  of  as  a  density- 

'  Pi+Pi+1  Pi+Pi+1 

weighted  face  velocity,  and  Pi+1/2  =  is  the  cell  face  density.  Note  that 
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we  currently  use  {pu)i^R  =  {pu)i  and  {pu)i+i^L  =  {.pu)i+ii  although  higher  order 
approximations  could  be  used.  Using  this  discretization  on  equation  (10)  yields 


1.2/^T 


I  +  p^{cT^t^G 


n+l 


G 
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where  G  is  our  gradient  discretization  and  — is  our  divergence  discretization 
and  the  hat  variables  are  defined  as  above.  This  is  solved  to  obtain  at 
cell  centers. 


It  is  interesting  to  note  that  this  derivation  does  not  require  an  ideal  gas 
assumption,  and  hence  should  be  general  enough  to  work  with  any  equation 
of  state  (even  multi-species  flow  [2]). 


2.2  Updating  Momentum  and  Energy 


To  obtain  the  correct  shock  speeds  we  need  the  pressure  at  cell  faces  for 
equations  (2)  and  (3),  and  the  velocity  at  cell  faces  for  equation  (3).  Applying 
conservation  of  momentum  to  the  control  volumes  Gi^r  and  (see  figure  1) 

gives 

Dui^R/Dt  =  {pi  -  Pi+112) / {/^xpi^R/2) 

and 

Dui+i^i/Dt  =  (pi+i/2  -Pi+i)/(Aa;pi+i,i/2). 

The  constraint  that  the  interface  remain  in  contact  implies  that  Dui^R/Dt  = 
Dui+i^l/ Dt,  which  can  be  used  with  the  aforementioned  equations  to  solve  for 
the  pressure  at  the  flux  location  Aj_|_i/2  as 


Pi+l/2 


Pi+lPi  +  PiPi+l 
Pi+1  +  Pi 


(15) 


For  solid  wall  boundaries,  we  reflect  the  pressure  and  density  values  as  usual, 
and  then  use  equation  (15).  The  cell  face  velocity  is  computed  via  equa¬ 
tion  (13),  and  pj+i/2hi+i/2  is  used  in  equation  (3). 


3  Time  Step  Restriction 


The  eigenvalues  of  the  Jacobian  of  the  advection  part  of  the  flux  are  all  u. 
Since  we  solve  the  acoustic  component  implicitly,  we  no  longer  have  a  severe 
time  step  restriction  determined  by  the  speed  of  sound  c,  and  all  that  remains 
is  to  find  an  estimate  for  the  maximum  value  of  |u|  throughout  the  time  step. 
Simply  using  u'^  is  not  enough,  since  e.g.  Sod  shock  tube  starts  out  with  an 
initial  velocity  identically  zero  and  thus  would  imply  an  infinite  At.  To 
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alleviate  this,  we  add  a  term  that  estimates  the  change  in  velocity  over  a  time 
step  similar  to  what  was  done  in  [8].  This  requires  consideration  of  which 
we  include,  in  our  estimate  of  the  velocity  at  the  end  of  the  time  step  to  get 

^  and  the  CFL  condition  becomes 
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which  is  quadratic  in  At  with  solutions 


(16) 
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As  the  lower  limit  is  always  non  positive  and  At  >  0,  we  only  need  enforce 
the  upper  bound.  As  px  — >  0,  both  the  numerator  and  denominator  vanish 
and  thus  we  obtain  a  more  convenient  time  step  restriction  by  replacing  the 
2nd  equation  (16)  with  this  upper  bound  to  obtain 


At 
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and  our  hnal  CFL  condition  becomes 
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Note  that  this  is  not  linear  in  Ax,  but  as  Ax  — 0  we  obtain  a  more  typical 
CFL  condition  At  <  — • 


4  Modified  ENO  Scheme 


When  using  a  traditional  ENO  methods  for  the  advection  part  of  our  equa¬ 
tions  (as  in  [11]),  we  obtained  excessive  spurious  oscillations.  This  seems  to 
be  related  to  our  dual  cell  center  and  MAC  grid  formulation,  thus  we  de¬ 
vice  a  new  ENO  scheme  which  better  utilizes  that  dual  formulation.  We  call 
this  Mach-ENO  or  MENO.  The  main  idea  is  to  replace  the  advection  velocity 
with  the  MAC  grid  value  defined  at  the  flux  in  question,  i.e.  it.  The  lowest 
level  of  the  divided  difference  table  is  typically  constructed  with  the  physical 
fluxes,  i.e.  pw,  pu^  and  Eu  for  Fi(U)  in  equation  (1).  A  dissipation  term  is 
added  for  the  local  and  global  Lax-Friedrichs  versions.  Consider  constructing 
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(a)  Standard  ENO-LLF  (b)  MENO-LLF 

Fig.  2.  Sod  shock  tube  problem  at  t  =  .15s.  Left:  Standard  ENO-LLF  using  401  grid 
points  (green)  and  1601  grid  points  (red).  Right:  The  base  1601  grid  points  solution 
is  the  same  as  in  the  left  figure,  but  the  coarse  grid  calculation  (with  401  grid  points) 
is  done  with  the  new  MENO  scheme.  Both  simulations  were  done  with  explicit  time 
stepping  and  a  full  characteristic  decomposition  in  order  to  demonstrate  that  the 
new  ENO  schemes  performs  similar  to  the  old  one  when  one  is  not  using  our  new 
implicit  discretization  of  the  pressure. 

an  ENO  approximation  for  the  flux  at  Xj+i/2-  Locally,  we  would  use  a  divided 
difference  table  with  base  values  corresponding  to  the  physical  fluxes  plus  or 
minus  the  appropriate  dissipation.  Our  modihcation  is  to  replace  PjUj,  PjUj"^, 
and  EjUj  with  pjMj+1/2,  PjUjUi+1/2,  and  EjUi+1/2  leaving  the  dissipation  terms 
unaltered.  Note  that  Ui+1/2  is  hxed  throughout  the  divided  difference  table 
similar  to  the  way  one  Exes  the  dissipation  coefficient. 

Figure  2  compares  our  new  MENO  scheme  to  the  standard  scheme  from  [11] 
for  standard  Sod  Shock  tube.  For  this  problem  the  results  are  fairly  similar, 
but  for  other  test  cases  the  MENO  scheme  performed  much  better  and  in 
fact  the  standard  ENO  scheme  was  not  successful  in  producing  any  solution 
whatsoever  for  hgure  10  in  our  examples  section. 


5  Numerical  Results 


We  use  third  order  TVD  Runge-Kutta  [10]  for  all  our  examples.  We  use  two 
variations  of  the  third  order  TVD  Runge-Kutta  scheme  in  all  of  our  examples. 
The  hrst  is  to  perform  Runge-Kutta  on  just  the  advection  part,  Fi(U),  with 
only  one  hnal  implicit  solve  for  F2(U).  The  second  variation  is  to  carry  out 
both  Fi(U)  and  F2(U)  for  each  Runge-Kutta  stage  although  this  has  three 
times  the  computational  cost  as  far  as  the  implicit  solution  of  F2(U)  is  con- 


cerned.  However  better  numerical  results  are  obtained  using  the  implicit  solve 
within  each  stage  of  the  Runge-Kutta  cycle  (see  hgure  3),  and  thus  that  is  the 
scheme  employed  throughout  the  rest  of  the  example  section,  unless  otherwise 
mentioned. 


(a)  One  implicit  solve  (b)  Three  implicit  solves 


(c)  One  implicit  solve  (d)  Three  implicit  solves 

Fig.  3.  Numerical  results  comparing  placing  the  implicit  solve  either  inside  each 
Runge-Kutta  stage  (b  and  d)  or  once  after  a  full  three  stage  Runge-Kutta  cycle 
(a  and  c).  The  top  two  figures  show  the  results  for  a  Sod  shock  tube  problem  at 
t  =  .15s,  the  bottom  two  figures  show  the  results  for  a  strong  shock  tube  problem 
at  t  =  2.5  X  10“®s.  Density  is  shown  in  all  figures.  Note  the  spurious  overshoots 
when  the  implicit  solve  is  not  included  in  the  Runge-Kutta  cycle  (left  two  figures). 
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5.1  One  dimensional  Validation 


For  the  one  dimensional  tests,  we  use  a  computational  domain  of  [0, 1],  401 
grid  points,  and  also  plot  a  baseline  solution  using  1601  grid  points  in  the 
standard  fully  explicit  ENO  method  as  in  [11],  A  second  order  ENO  was  used 
along  with  the  CEL  number  of  .5. 


5.1.1  Sod  Shook  Tube 


Our  Erst  test  case  is  a  standard  Sod  shock  tube  with  initial  conditions  of 


(p(a;,0),M(a;,0),p(a;,0)) 


(1,0,1)  if  a;  <.5, 

(.125,0,.!)  if  a;  >.5. 


Our  results  are  shown  in  Figure  4,  which  indicate  well  resolved  shock,  rar¬ 
efaction  and  contact  solutions.  Since  our  method  is  conservative,  we  get  the 
correct  shock  speeds. 


5.1.2  Lax’s  Shook  Tube 


Lax’s  shock  tube  is  similar  in  nature  to  Sod  shock  tube,  except  that  the  initial 
condition  has  a  discontinuity  in  the  velocity. 


(p(a;,0),M(a;,0),p(a;,0)) 


(.445,  .698,  3.528)  if  a;  <  .5, 
(.5,0,  .571)  if  a;  >.5. 


Our  results  are  shown  in  Figure  5. 


5.1.3  Strong  Shock  Tube 

The  Strong  shock  tube  problem  poses  initial  conditions  that  generates  a  su¬ 
personic  shock. 


(p(a;,0),M(a;,0),p(a;,0)) 


(1,0, 10^°)  if  a;  <.5, 
(.125,0,.!)  if  a;  >.5. 


Our  results  are  shown  in  Figure  6.  However  the  scheme  admits  some  oscilla¬ 
tions  near  the  rarefaction  wave.  Note  that  the  main  advantage  of  the  proposed 
method  is  to  take  target  time  steps  irrespective  of  the  sound  speed  values,  one 
could  use  the  usual  ENO  scheme  for  high  Mach  number  flows  (or  high  Mach 
number  regions  of  the  flow  -  if  asynchronous  timestepping  is  used). 
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(a)  Density  (b)  Velocity 


(c)  Pressure  (d)  Internal  Energy 

Fig.  4.  Numerical  Results  of  the  Sod  shock  tube  problem  at  t  =  .15s.  The  explicit 
baseline  solution  is  plotted  in  red,  and  the  solution  from  our  method  is  plotted  in 
dotted  green. 


5.1.4  Mach  3  Shock  Test 

The  initial  conditions  for  the  Mach  3  shock  test  are: 


(p(a;,0),M(x,0),p(a;,0)) 


(3.857,  .92, 1.333)  if  x  <  .5, 
(1,3.55,1)  if  a;  >.5. 


Our  results  are  shown  in  Figure  7.  As  above  we  do  note  some  oscillations  near 
the  rarefaction  wave. 


11 


(a)  Density 


(b)  Velocity 


(c)  Pressure  (d)  Internal  Energy 

Fig.  5.  Numerical  Results  of  the  Lax’s  shock  tube  problem  at  t  =  .12s.  The  explicit 
baseline  solution  is  plotted  in  red,  and  the  solution  from  our  method  is  plotted  in 
dotted  green. 


5.1.5  High  mach  flow  test 

The  initial  conditions  for  the  High  mach  flow  test  are: 


(p(a;,0),M(a;,0),p(a;,0)) 


(10,2000,500)  if  a:  <.5, 
(20,0,500)  ifx>.5. 


As  noted  in  [7]  the  Mach  number  in  this  test  can  reach  as  high  as  240.  Our 
results  are  shown  in  Figure  8. 
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(a)  Density 


(b)  Velocity 


(c)  Pressure  (d)  Internal  Energy 

Fig.  6.  Numerical  Results  of  the  strong  shock  tube  problem  at  t  =  2.5  x  10“®s. 
The  explicit  baseline  solution  is  plotted  in  red,  and  the  solution  from  our  method 
is  plotted  in  dotted  green. 

5.1.6  Interaction  of  blast  waves 

Here  we  present  a  test  of  two  interacting  blast  waves.  This  problem  was  intro¬ 
duced  by  [12]  and  involves  multiple  strong  shock  waves.  The  initial  conditions 
for  the  test  are: 


[(1,0,103)  if0<a:<.l, 

{p{x,  0),u{x,  0),p{x,  0))  =  <  (1,  0, 10“^)  if  .1  <  X  <  .9, 

[(1,0,102)  if.9<x<l. 


We  also  have  solid  wall  boundary  conditions  at  x  =  0  and  x  =  1.  Our  results 
are  shown  in  Figure  9  which  shows  that  we  achieve  very  accurate  results. 
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(a)  Density  (b)  Velocity 


(c)  Pressure  (d)  Internal  Energy 

Fig.  7.  Numerical  Results  of  the  Mach  3  shock  tube  problem  at  t  =  .09s.  The  explicit 
baseline  solution  is  plotted  in  red,  and  the  solution  from  our  method  is  plotted  in 
dotted  green. 

5.1.7  Two  Symmetric  Rarefaction  Waves 

In  this  test  there  are  two  rarefaction  waves  going  in  opposite  directions  from 
the  center  of  the  domain.  This  causes  very  low  density  regions  near  the  center 
of  the  domain.  The  initial  conditions  for  the  test  are: 


(p(x,0),M(x,0),p(a:,0)) 


(1,-2,  .4)  ifx<.5, 
(1,  2,  .4)  if  X  >  .5. 


Our  results  are  shown  in  Figure  10.  Our  results  are  comparable  to  that  of  [7] 
and  [13].  Note  that  there  is  an  unphysical  pulse  in  the  internal  energy  field 
near  the  low  pressure  region,  caused  by  overheating  (see  e.g.  [3]). 


14 


(a)  Density 
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(c)  Pressure  (d)  Internal  Energy 

Fig.  8.  Numerical  Results  of  the  High  Mach  shock  tube  problem  at  t  =  1.75  x 
The  explicit  baseline  solution  is  plotted  in  red,  and  the  solution  from  our  method 
is  plotted  in  dotted  green. 

5.1.8  Smooth  Flow  Test  (Mach  Zero  Limit) 


The  initial  conditions  for  the  zero  mach  limit  test  are  given  by: 


u{x,  0)  =  0 
pix,0)  =po  +  epi{x) 

Pi{x)  =  60cos(27rx)  +  100sin(47rx) 


p{x,0) 


(^)  ’ 


Where  po  =  1,  Po  =  10®,  and  e  =  10^.  Since  the  flow  is  smooth  and  there  are 
no  shocks  in  this  test,  we  have  used  a  single  implicit  solve  per  time  step.  This 
test  is  dominated  by  acoustic  waves  (as  observed  in  [7]).  We  can  take  time 
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(a)  Density  (b)  Velocity 


(c)  Pressure  (d)  Internal  Energy 

Fig.  9.  Numerical  Results  of  the  interacting  blasts  shock  tube  problem  at  t  =  .038s. 
The  explicit  baseline  solution  is  plotted  in  red,  and  the  solution  from  our  method 
is  plotted  in  dotted  green. 

steps  as  large  as  is  permitted  by  our  CFL  condition  in  equation  (17).  This 
permits  time  steps  three  orders  of  magnitude  greater  than  those  permitted 
by  sound-speed  based  CFL.  However,  as  with  all  implicit  schemes,  taking  too 
large  a  time  step  can  lead  to  inaccurate  results.  Thus,  in  order  to  get  sufficient 
accuracy,  we  clamp  our  time  step  to  be  a  hxed  multiple  of  the  explicit  time 
step.  In  hgure  11  we  use  3  times  the  explicit  time  step  and  show  convergence 
via  grid  resolution.  In  a  second  suit  of  tests  we  show  that  we  can  increase 
the  grid  resolution  without  the  need  for  rehning  the  time  step.  The  timing 
results  for  this  experiment  are  available  in  Table  13  where  At  remains  hxed 
as  the  grid  resolution  goes  up  as  high  as  320,  000  gird  cells.  At  that  point  the 
effective  sound  speed  CFL  is  300.  Numerical  results  are  plotted  in  hgure  12 
and  Table  13  summarizes  the  results.  In  particular  whereas  the  newly  proposed 
implicit  method  permits  a  hxed  time  step  all  the  way  up  to  320,  000  grid  points 
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(a)  Density 


2 
1.5 
1 

0.5 
0 

-0.5 
-1 
-1.5 
-2 

0  0.2  0.4  0.6  0.8  1 


(b)  Velocity 
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(c)  Pressure  (d)  Internal  Energy 

Fig.  10.  Numerical  Results  of  the  symmetric  rarefaction  shock  tube  problem  at 
t  =  .15s.  The  explicit  baseline  solution  is  plotted  in  red,  and  the  solution  from  our 
method  is  plotted  in  dotted  green. 

allowing  the  wall  clock  simulation  time  to  scale  approximately  linear  to  the 
size  of  the  problem,  the  explicit  simulation  time  grows  quadratically  even 
becoming  impractical  at  320,  000  grid  points. 


5.2  Flow  Past  a  Step  Test 


Our  hrst  two  dimensional  experiment  is  similar  to  the  one  described  in  [3]. 
We  assume  an  ideal  gas  with  7  =  1.4.  The  test  domain  is  3  units  long  and  1 
unit  wide,  with  a  .2  unit  high  step  which  is  located  .6  units  from  the  left  hand 
side  of  the  tunnel.  The  initial  conditions  are  p  =  1.4,  p  =  1  and  u  =  3  and 
n  =  0  everywhere  in  the  domain.  We  apply  an  inflow  boundary  condition  on 
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Frame  15 


Fig.  11.  Numerical  results  comparing  the  pressure  in  smooth  flow  test  at  200,  400, 
800,  1600,  and  3200  grid  cells  with  an  effective  sound  speed  based  CFL  number  3 
at  t  =  1.5  X  10“®s.  The  red  curve  is  the  explicit  simulation  run  at  3200  grid  cells 
with  a  CFL  number  .5. 

the  left  hand  side  of  the  domain,  and  an  ontflow  bonndary  condition  on  the 
right  hand  side  of  the  domain.  A  reflective  solid  wall  bonndary  condition  is 
applied  for  the  top  and  bottom  bonndaries  of  the  domain.  We  show  nnmerical 
resnlts  at  f  =  4s  on  a  grid  of  resolntion  120x40  in  hgnre  14. 


5.3  Circular  Shock  Test 


The  circnlar  shock  test  has  an  initial  condition  prescribed  as 


{p,u,v,p) 


(1,0, 0,1)  ifr<.4 
(.125,0,0,.!)  ifr>.4 
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Frame  15 


Fig.  12.  Numerical  results  showing  pressure  in  smooth  flow  test  at  3200,  32000  and 
320000  grid  cells.  We  used  an  effective  sound  speed  based  CFL  number  of  3,  30  and 
300  respectively  at  t  =  1.5  x  10“^s.  Since  At  stays  constant,  the  solution  remains 
relatively  unchanged  even  as  we  get  huge  time  step  gains. 
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sound  speed 

CFL 

At 
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Wall  clock  time 

(Explicit) 
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511.67s 

32000 
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5.01-08 

810.03s 

60498.49s 

320000 

300 

5.01e-08 

9976.58s 

Impractical 

Fig.  13.  Timing  results  for  smooth  flow  test,  with  At  approximately  constant.  The 
wall  clock  times  are  shown  for  simulations  till  t  =  5  x  10“^s. 


where  r  =  +  y'^.  Numerical  results  are  shown  in  hgure  15.  The  same 

test  was  shown  in  [14].  Our  results  indicate  well  resolved  shock  and  contact 
solutions  along  with  correct  speed  shock  calculations. 
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Fig.  14.  Numerical  results  showing  the  Schlerian  plots  of  density  for  the  flow  past  a 
step  test  on  a  grid  of  size  120x40  at  t  =  4s. 


(a)  Density  (b)  Pressure 

Fig.  15.  Numerical  results  for  the  circular  shock  test  on  a  grid  of  size  100x100  at 
t  =  .25s. 


6  Conclusions  and  Future  Work 


We  have  presented  a  method  for  alleviating  the  stringent  CFL  condition  im¬ 
posed  by  the  sonnd  speed  in  highly  non-linear  compressible  flow  simulations. 
A  fractional  step  procedure  combined  with  the  pressure  evolution  equation 
is  used.  The  method  works  for  arbitrary  equations  of  state;  and  in  the  limit 
as  the  sound  speed  goes  to  inhnity,  it  yields  the  Poisson  equation  for  incom¬ 
pressible  flow.  We  also  presented  a  Mach-ENO  or  MENO  scheme  which  better 
utilizes  a  dual  cell  center  and  MAC  grid  formulation.  The  numerical  exper¬ 
iments  on  various  benchmark  problems  for  one  and  two  dimensions  indicate 
that  our  semi-implicit  method  obtains  well  resolved  shock,  rarefaction  and 
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contact  solutions.  Since  our  method  is  conservative,  we  also  obtain  correct 
shock  speeds.  The  smooth  flow  example  illustrates  the  ability  of  our  method 
to  take  significantly  large  time  steps  for  low  Mach  number  flows  as  compared 
to  explicit  methods.  In  future  work  we  plan  to  extend  our  approach  to  handle 
two-way  coupling  between  compressible  and  incompressible  flows,  as  well  as 
fully  implicit  solid-fluid  coupling. 
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A  Boundary  Conditions 


Figure  14  requires  the  handling  of  inflow  and  outflow  boundary  conditions. 
We  define  Uout  to  be  the  outgoing  state  and  Uin  to  be  the  ingoing  state.  The 
outgoing  state,  Uout,  is  obtained  by  simple  extrapolation  whereas  the  ingoing 
state,  Uin,  is  obtained  by  attenuating  Uout  towards  specified  far-fleld  values. 
After  defining  Uout  via  extrapolation,  we  average  the  primitive  variables  to 
cell  flux  on  the  boundary  of  the  domain,  and  use  those  values  to  compute  a 
characteristic  decomposition.  If  the  characteristic  held  indicates  ingoing 
information,  then  when  applying  the  ENO  scheme  in  this  characteristic  held 
we  use  Uin  for  the  ghost  node  values.  Otherwise  Uout  is  used.  Note  for  higher 
order  schemes  boundary  values  will  be  needed  for  fluxes  on  the  interior  of  the 
domain  as  well,  and  we  choose  the  ghost  nodes  (as  Uin  or  Uout)  in  the  same 
fashion. 

Our  ingoing  state,  Uin,  is  calculated  as  follows.  Our  ingoing  state,  Uin,  is  ob¬ 
tained  by  attenuating  the  extrapolated  state,  Uout,  towards  a  given  far  held 
state,  Ufar-  This  is  accomplished  by  multiplying  Uout  with  each  of  the  left 
eigenvectors,  attenuating  if  the  eigenvalue  in  that  characteristic  held  indicates 
an  ingoing  wave,  and  then  multiplying  by  the  right  eigenvector.  Defining  the 
scalar  characteristic  information  in  each  held  as  =  L^Uout,  we  would  atten¬ 
uate  towards  using  the  analytic  solution  of  the  ODE 

d^/dt  =  K{^-^for) 
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for  time  step  At  using  initial  data  of  .^  =  ^out-  We  used  an  attenuation  coeffi¬ 
cient  of  =  —  .5  in  our  examples. 
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